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ABSTRACT 


We explore a possible scenario of the explosion as a result of core collapses of rotating massive 
stars that leave a black hole by performing a radiation-viscous-hydrodynamics simulation in numerical 
relativity. We take moderately and rapidly rotating compact pre-collapse stellar models derived in 
stellar evolution calculations as the initial conditions. We find that the viscous heating in the disk 
formed around the central black hole powers an outflow. For rapidly rotating models, the explosion 
energy is Z 3 x 105! erg, which is comparable to or larger than that of typical stripped-envelope 
supernovae, indicating that a fraction of such supernovae may be explosions powered by black-hole 
accretion disks. The explosion energy is still increasing at the end of the simulations with a rate of 
> 109? erg/s, and thus, it may reach ~ 10° erg. The nucleosynthesis calculation shows that the mass of 
56Ni amounts to > 0.1 Mo, which, together with the high explosion energy, satisfies the required amount 
for broad-lined type Ic supernovae. The moderately rotating models predict small ejecta mass of order 
0.1 M and explosion energy of < 10°! erg. Due to the small ejecta mass, these models may predict a 
short-timescale transient with the rise time 3-5 d. It can lead to a bright (~ 10“4 erg/s) transient like 
superluminous supernovae in the presence of dense massive circum-stellar medium. Irrespective of the 
models, the lowest value of the electron fraction of the ejecta is = 0.4, and thus, the synthesis of heavy 
r-process elements is not found in our calculation. 


Keywords: stars: neutron; general-hydrodynamics-neutrinos-relativistic processes 


1. INTRODUCTION with no successful shock revival due to the strong ram 
pressure by the matter infall (e.g., O'Connor & Ott 
2011). Even in this case, there is a possibility of the ex- 
plosion if a massive disk is formed around the black hole 
due to the rotation of the progenitor star. In the disk, 
the magnetorotational instability (MRI) could amplify 
the magnetic field, developing a turbulent state inside 
the disk (e.g., Balbus & Hawley 1991). The turbulent 
motion then induces an effective viscosity in the disk, 


At the final evolution stage of massive stars, their iron 
cores become gravitationally unstable and collapse. Af- 
ter the core bounce due to the formation of a proto- 
neutron star, a shock wave is formed on its surface and 
propagates outward. The shock wave then stalls pri- 
marily because of the photodissociation of heavy (iron- 
group) nuclei in the infalling matter swept by the shock. 


In the standard neutrino-driven delayed-explosion sce- : t 
nario, the shock is revived by the heating of neutrinos which governs the evolution of the disk through the an- 


emitted by the proto-neutron star (e.g., Janka et al. gular momentum transport and heating. The viscous 
2012). heating rate is estimated by 


If the collapsing star has a very compact core, the Lys ~ VO? Mai 


proto-neutron star is likely to collapse into a black hole 3 
& 7 x 109. erg/s ( “8 = 
0.03 / \ 109 cm/s 
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where Mgg is the mass of the black hole, Maia; and 
Tdisk are the mass and typical radius of the disk, c, is 
the sound speed, Q is the local angular velocity, and we 
assumed a Kerlerian rotation and the Shakura-Sunyaev- 
type alpha-viscosity model for the kinetic viscous coef- 
ficient v = ayiscsH (Shakura & Sunyaev 1973) with the 
disk scale height H = c,/Q. oi is the so-called al- 
pha parameter, which is likely to be of order 107? in 
the presence of MRI turbulence (e.g., Balbus & Hawley 
1998; Hawley et al. 2013; Suzuki & Inutsuka 2014; Shi 
et al. 2016; Kiuchi et al. 2018). Accretion disks formed 
around black holes in the collapsar scenario (MacFadyen 
& Woosley 1999) often become a neutrino-dominated 
accretion disk (e.g., Kohri & Mineshige 2002; Di Mat- 
teo et al. 2002) in which internal energy generated by 
the viscous heating is released primarily by the neutrino 
emission. However, the latest studies for the systems of 
a black hole and à compact accretion disk have shown 
that in the late stage of the viscous evolution, the neu- 
trino cooling late drops due to the viscous angular mo- 
mentum transport and subsequent expansion of the disk 
(Fernández & Metzger 2013; Just et al. 2015; Fujibayashi 
et al. 2020a,b; Just et al. 2022a). In such a stage, the 
viscous heating can be used for launching a strong out- 
flow from the disk in a viscous timescale estimated by 


r2, Bs C E 
Luis "T disk x As vis s 
v 0.03 10° cm/s 


Mpu VU ras XS 
x . (2) 
10M 107 cm 


The order of magnitude of the explosion energy gener- 
ated by the viscous heating is estimated by ~ Lyistvis, 
which is comparable to or even larger than that of typ- 
ical supernovae (~ 10°! erg) for plausible values of ois, 
Mpu, Maisk, fdisk, and cs. This motivates us to explore 
a scenario of the explosion from a massive accretion disk 
around a spinning black hole formed during the rotating 
stellar core collapse. 

'The sub-relativistic outflow from the disk is of impor- 
tance for several aspects. First, it can be an impor- 
tant energy source to power a supernova-like explosion 
associated with long-duration gamma-ray bursts (e.g., 
MacFadyen & Woosley 1999, Pruet et al. 2003, Na- 
gataki et al. 2007, Surman et al. 2006, and Hayakawa 
& Maeda 2018), for which the promising central engine 
of the gamma-ray bursts is likely to be a spinning black 
hole penetrated by a strong magnetic field. However, 
additionally, we need the supernova component. Eisen- 
berg et al. (2022) suggested, based on their simulations, 
that the observationally inferred velocity distribution of 
the supernova-component is not likely reproduced only 
by the relativistic jet. This indicates that there has to 


be another energy source to drive the supernova com- 
ponent in addition to the relativistic jet accounting for 
the gamma-ray burst. Kohri et al. (2005) applied a 
disk explosion scenario to normal supernovae by analyti- 
cally solving a stationary neutrino-cooled accretion-disk 
model and indicated that an energetic outflow could be 
driven from the collapse of rotating stars when the ac- 
cretion flow is advection-dominated.! 

Second, it has been speculated that the matter in the 
disk outflow could be neutron-rich, and thus, the out- 
flow may be a site for r-process nucleosynthesis (Surman 
et al. 2006; Pruet et al. 2003; Kohri et al. 2005). Siegel 
et al. (2019) suggested, based on their magnetohydrody- 
namics simulations with an approximate neutrino treat- 
ment in a fixed black-hole spacetime, that neutron-rich 
matter may be ejected from the disk cooled by neutrinos 
and the heavy nuclei up to third peak of the r-process el- 
ements may be synthesized. In a similar setup but with 
Monte-Carlo neutrino transfer, however, Miller et al. 
(2020) pointed out that the electron fraction (Ye) of the 
ejecta is higher than 0.3, and thus, nuclei only up to sec- 
ond peak of the r-process elements are synthesized. Just 
et al. (2022b) performed viscous hydrodynamics simula- 
tions in Newtonian gravity with general relativistic cor- 
rections incorporating moment-based neutrino radiation 
transfer for the collapse of a rotating massive star and 
showed that the outflow from the disk formed around 
the black hole has the electron fraction higher than 0.4. 
Therefore, the speculation in this field has not converged 
yet. Moreover, no fully general relativistic work, which 
self-consistently takes into account the self gravity of the 
collapsing star and the formed black hole in general rel- 
ativistic manner, has been carried out. Obviously more 
detailed studies are required. 

Third, recent high-cadence transient surveys have 
shown that there is a variety of optical transients that 
are not canonical supernovae. Those with timescales 
of a few days, which are much shorter than that of 
normal supernovae (> 10d), are such examples (e.g., 
Drout et al. 2014, Prentice et al. 2018, and Tampo et al. 
2020). Despite intensive photometric and spectral obser- 
vations, the progenitors of the transients different from 
the canonical supernovae are still not clear. There are 
several scenarios in which a collapse of a massive star 
leading to black-hole formation plays a central role (e.g., 
Margutti et al. 2019, Perley et al. 2019). However, the 


1 Note that it is still possible to achieve successful neutrino-driven 
explosion in proto-neutron star phases even for collapses for com- 
pact progenitor stellar cores (see a series of work Obergaulinger 
& Aloy 2020; Aloy & Obergaulinger 2021; Obergaulinger & Aloy 
2021, 2022 and also Fujibayashi et al. 2021). 
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previous works are limited only to the ones based on sim- 
plified models (see Piran et al. 2019 and Gottlieb et al. 
2022b for a recent simulation-based model). Thus, it is 
important to provide predictions based on reliable nu- 
merical simulations for interpreting the observation and 
confirming the origins of the mysterious transients. 

Motivated by these current situations, in this pa- 
per, we explore the long-term evolution of the collapse 
of rotating massive stars by fully general relativistic 
radiation-viscous-hydrodynamics simulations with an 
approximate neutrino transfer. The numerical simula- 
tions are performed employing two types of the initial 
condition. The first one is based on a non-rotating pre- 
collapse stellar model evolved from a 20M, helium star 
(Takahashi et al. 2018), for which we add a stellar rota- 
tion. The other is based on rapidly rotating pre-collapse 
stellar models calculated by Aguilera-Dena et al. (2020). 
For these pre-collapse stellar models, we explore how the 
disk is formed and evolved around the black hole dur- 
ing rotating stellar collapse, leading to a supernova-like 
stellar explosion. 

'This paper is organized as follows. In 8 2, our method 
of the simulations is briefly described. We also intro- 
duce the pre-collapse stellar models which we employ 
from the stellar evolution calculations. Then, in § 3, the 
results of our numerical-relativity simulations and the 
nucleosynthesis calculations are presented. We discuss 
the possible optical transients and implications to broad- 
lined type Ic supernovae and gamma-ray bursts based on 
our results in 8 4. We also discuss the possible produc- 
tion of light r-process nuclei and effects on the optical 
transient. 8 5 is devoted to a summary. Throughout this 
paper, G, c, and kg denote the gravitational constant, 
speed of light, and Boltzmann's constant, respectively. 


2. METHOD 
2.1. Numerical code 


Numerical-relativity simulations are performed with 
our latest axisymmetric neutrino-radiation viscous- 
hydrodynamics code. The detail of the code is described 
in Fujibayashi et al. (2017, 2020c). In this code, Ein- 
stein's equation is solved in the original version of the 
Baumgarte-Shapiro-Shibata-Nakamura formalism (Shi- 
bata & Nakamura 1995; Baumgarte & Shapiro 1998) 
with a constraint propagation prescription to make 
the constraint violation to propagate outward (Hilditch 
et al. 2013). A dynamical gauge condition described in 
Fujibayashi et al. (2017) is employed. To impose the 
axisymmetry for the geometrical variables, the so-called 
cartoon method (Alcubierre et al. 2001; Shibata 2000) 
with the fourth-order Lagrange interpolation is imple- 
mented. 


The neutrino radiation transfer equations are approx- 
imately solved using a leakage scheme together with the 
truncated moment formalism (Fujibayashi et al. 2017; 
see also Sekiguchi 2010). In this formulation, the neu- 
trino field is split into two components; trapped and 
free-streaming neutrinos. The trapped neutrinos are as- 
sumed to be tightly coupled with the fluid and have the 
same local temperature and velocity as those of the fluid. 
This component is treated as a part of the fluid and con- 
tributes to the internal energy and pressure. It becomes 
the free-streaming component with the generation rate 
controlled by the local diffusion rate of neutrinos. The 
free-streaming neutrinos are assumed to obey radiation 
transfer equations, which are solved by a truncated mo- 
ment formalism with the M1 closure relation (Thorne 
1981; Shibata et al. 2011). Following our previous work 
(Fujibayashi et al. 2020a,b,c), we solve the equations for 
the frequency-integrated energy and momentum density 
for the three neutrino radiation fields (electron, electron 
anti-, and other neutrinos). 

The viscous hydrodynamics equations are solved using 
the formulation described in Shibata et al. (2017), in 
which the energy-momentum tensor is written as 

Tj, = phu,u, + Dën, + vphro,, (3) 
where p :— Mum is the rest-mass density, h = c? + € + 
P/p is the specific enthalpy with the specific internal 
energy € and pressure P, and v is the kinematic viscous 
coefficient. The viscous tensor 77, is a symmetric tensor 
that satisfies (Israel & Stewart 1979) zl up = 0, and it 
is determined by the following equation 


Lon = E Se Gul, (4) 
where Lu denotes the Lie derivative with respective to 


u” and C is a coefficient, for which we set CT = O(10ys). 
Assuming the form of the shear tensor as 


gu = V ptv + Vyuy = Lu(Guv + upt), (5) 


we obtain the evolution equation for Tuy :— fu —G(gu» E 
UpUy) as 


af zs e (6) 
We only need to solve the spatial part of zu, because 
of the presence of the condition gu" — 0. The spatial 


part obeys the following evolution equation (in Carte- 
sian coordinates): 


O (pu! Griz) + Ox (pu^ V/—gri;) 
+ pu / —g(rixOjv^ + rj Ov") = -ev-g6Cr. (7) 
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Table 1. Model description. The columns provide from left to right: model name, progenitor model name, mass of the 
progenitor star, angular velocity profile, equation of state, compactness just prior to the collapse, innermost grid spacing, the 
values of ô and N, and the location of the outer boundaries along each axis. The last column shows a model for the black 


hole-disk system. 


Model Progenitor Progenitor star Q profile EOS €25 Azo (m) A N L (cm) 
T20 T20* Mue = 20M Eq.(10)  SFHo 0.74 150 0.0 991 1.0x 107? 
ADO9xi ` ADO9^ Mzams = 9Mo Original x1 DD2 0.68 175 0.01 975 1.0x 10!? 
AD20xi  AD20° Mzams = 20Mo Original x1 DD2 0.66 175 0.01 975 1.0 x 10!? 
AD20x2 ` AD20^ Original x2 

T201 T20* lower-res. T20 SFHo 0.74 150 0.013 807 1.0 x 107? 
ADO9xil  ADO9P lower-res. ADO9x1 DD2 0.68 175 0.013 793 1.0 x 10!? 
BHdisk  — 10M BH-3Mo disk = DD2 - 220 0.01 801 3.6 x 10? 


Takahashi et al. (2018), PAguilera-Dena et al. (2020) 


'The effective viscosity in the disk is believed to arise as 
a result of the turbulence induced by magnetohydrodyn- 
mical instabilities such as MRI (Balbus & Hawley 1991; 
Balbus & Hawley 1998) and Kelvin-Helmholtz instabil- 
ity (e.g., Obergaulinger et al. 2010). Following Shakura 
& Sunyaev (1973), we define the viscous coefficient by 


Hz CsÉtur; (8) 


where bur is the mixing length scale (or the largest eddy 
size) in the turbulence. In the alpha disk model, tur is 
written as avi H (Shakura & Sunyaev 1973). In this 
study, we assume 


2G Men Mpu 
Liur = 0.03 x — 7x 440 : 9 
5 5 c? Z6 ( ) 


where the black-hole mass, Mpy, continuously increases 
with time due to the matter accretion. Thus, we em- 
ploy a time-varying form of £44,. Since the scale height 
H is larger than 2GMpy/c? for most parts of the disk 
around the black hole, Equation (9) implies that we as- 
sume Qyis < 0.03, i.e., a conservative value of ayis. Even 
for such a conservative value of ayis, we will find a signif- 
icant effect in 8 3. We note that the viscosity is incorpo- 
rated just prior to the formation of the disk around the 
black hole. That is, for the formation and early growth 
stages of the black hole, we do not consider the viscous 
effects. 

The grid structure is the same as in the 2D sim- 
ulations recently performed with the same code (Fu- 
jibayashi et al. 2022), in which the cylindrical coordi- 
nates (R,z) are employed. In the inner cylindrical re- 
gion of R < 100A, and z € 100Azp9, a uniform grid 
with the grid spacing of Azo is prepared, while in the 
outer region, a non-uniform grid with an increase rate 
of the grid spacing of 1 + ô is prepared. The value of ó, 
grid number N for each axis, and location of the outer 
boundaries along each axis (denoted by L) are listed in 
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Figure 1. Specific angular momentum (solid curves) and 
that for innermost stable circular orbit (dashed curves) along 
the equatorial direction as functions of enclosed mass. 


Table 1. We assume the plane symmetry with respect 
to the z = 0 plane (the equatorial plane). 


2.2. models 


We employ the final states of massive stars evolved in 
stellar evolution calculations as the initial conditions of 
our numerical-relativity simulations. One of the progen- 
itors is the one evolved from a helium star with the initial 
mass of My. = 20Mo. This is taken from Takahashi 
et al. (2018) (hereafter T20). For this progenitor, the 
compactness at the radius with its enclosed mass 2.5Mo 
(referred to as r2.5 4), €2.5 :— 2.5/(r2.5 Mọ /1000 km), is 
0.74. 'This implies that it is very compact and expected 
to form a black hole with no neutrino-driven explosion if 
it is non-rotating (see, e.g., O'Connor & Ott 2011). Be- 
cause the progenitor is non-rotating, we add an artificial 
rotation profile as 


1 
7 1/Qo St 1/Qrk=ro (r)' 


Q(r) (10) 
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where Qo = 0.1 rad/s and 


MELIOR 

Qu Luv) NI (1) 
is the angular velocity for which the fluid element settles 
into a Keplerian motion with the cylindrical radius of ro. 
Note that the shellular rotation is assumed and the ini- 
tial angular velocity is a function only of the spherical 
radius r. In this study, we set ro = 1000 km. This im- 
plies that the maximum radius at which the infalling 
matter is circularized is 1000 km. In this model, we in- 
tentionally prepare a rapidly rotating progenitor star to 
investigate the effect of earlier disk formation around a 
formed black hole on the subsequent evolution of the 
system. 

'The other progenitors are rotating stars with zero-age- 
main-sequence masses of 9 and 20Mo, which are taken 
from Aguilera-Dena et al. (2020) (hereafter ADO9 and 
AD20, respectively). Because of their nearly chemically 
homogenious evolution, their pre-collapse profiles have 
very massive core and are very compact (£.5 = 0.68 and 
0.66 for ADO9 and AD20, respectively). The original ro- 
tation profile is applied for models ADO9x1 and AD20x1. 
For model AD20x2, on the other hand, the original angu- 
lar velocity is doubled to investigate the dependence of 
the result on the total angular momentum of the system. 

Figure 1 shows the specific angular momentum along 
the equatorial direction as a function of enclosed mass 
for the models employed in this paper. The specific an- 
gular momentum at the innermost stable circular orbit 
for the black hole with the enclosed mass and angular 
momentum, jrsco, is also plotted in the dashed curves. 
For each progenitor model, the enclosed mass at which 
the solid and dashed curves cross approximately gives 
the black-hole mass at which the infalling matter starts 
forming an accretion disk around the black hole. This 
figure indicates that disks are likely to be formed when 
the black holes grow to 4.0, 4.7, 4.9, and 7.8M for mod- 
els T20, AD20x2, ADO9x1, and AD20x1, respectively. In 
8 3, we will confirm that this prediction is approximately 
correct. 

At the time of the disk formation, the mass infalling 
rates for models T20 and AD20x2 are higher than those 
for ADO9x1 and AD20x1. This makes a difference in the 
evolution process of the disk and the onset timing of the 
outflow from the disk as found in § 3.2 and 8 3.3. 

In several previous studies, the mass ejection and 
nucleosynthesis calculation were carried out by simply 
modeling collapsar remnants as a black hole-disk sys- 
tem (Siegel et al. 2019; Miller et al. 2020). The evolu- 
tion process of the black hole-disk system as a result of 
stellar collapse may be different from that of an isolated 


disk around a black hole. To clarify the possible differ- 
ences, we also perform a simulation for a system of a 
spinning black hole and a massive disk (model BHdisk) 
as in our previous studies (Fujibayashi et al. 2020a,b). 
The initial condition is constructed using a code of Shi- 
bata (2007). The masses of the black hole and disk are 
10 and 3Mo, respectively, and the dimensionless spin 
of the black hole is ~ 0.6. The inner and outer edges 
of the disk are chosen to be 4GMpy/c? ~ 59km and 
400GMgu/c? ~ 5900km, respectively. The constant 
entropy per baryon of s = Tkg is simply assumed. 

In this paper, we employ tabulated equations of state 
(EOSs) referred to as SFHo (Steiner et al. 2013) and 
DD2 (Banik et al. 2014). For models with T20 pro- 
genitor, we employed SFHo EOS. On the other hand, 
we employed DD2 EOS for ADO9 and AD20 progenitors. 
Since the maximum mass of spherical neutron stars in 
the SFHo EOS is smaller than that in the DD2 EOS, a 
black hole is formed slightly earlier (i.e., the black-hole 
mass is smaller at its formation) in the model with the 
SFHo EOS. However, besides this difference, the depen- 
dence of the results on the EOS is weak because the 
maximum density of the disk is always smaller than the 
nuclear saturation density. 

To evolve black holes in a good precision, the radius 
of the apparent horizon has to be resolved well (for the 
dependence of the black-hole mass and spin on the grid 
resolution, see, e.g., Fujibayashi et al. 2020a). Employ- 
ing a stiller EOS, which predicts a larger maximum 
mass of neutron star, is advantageous to numerically 
resolve black holes in a better accuracy (for a given com- 
putational resource), because the black-hole mass (i.e., 
the radius of the apparent horizon) at its formation are 
larger for stiffer EOSs. Thus, for the present study, the 
DD2 EOS would be a better choice to save computa- 
tional resources. 


2.3. Diagnostic 


We define the mass infall rate to the central region by 


Mian = v/—gpu" dsp, (12) 


r=Tin 


where ds; = r75;-dcos @d@ is the surface element of a 
sphere and rj, is chosen in the following manner. In the 
early phase, it is the largest value among the radii of the 
surface of a standing accretion shock formed as a result 
of the core bounce. After this shock disappears due to 
the matter infall from the outer region which induces 
the collapse of a proto-neutron star to a black hole, the 
maximum radius of the apparent horizon ray is used for 
rig. After the formation of the disk, we again choose the 
largest value among the radii of the surface of a standing 
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'Table 2. Main results. The columns provide from left to right: model name, post-bounce time of black hole formation and 
explosion, ejecta mass, explosion energy, mass of the matter with maximum temperature higher than 5 GK (1GK — 10? K), 
mass of °°Ni, estimated peak bolometric luminosity, and rise time of bolometric light curve. For Mej and Eexp, there are two 
columns with different extraction radius. The values in parentheses for columns of Mej and Eexp are, respectively, the mass 
and binding energy of stellar matter above the extraction radius. Those for columns of Lpeak and trise are values for the case in 
which a half of the mass and binding energy above the extraction radius are considered. 


Model tBu — texp Mej Eexp Msscax Myi Lpeak trise 

Il (s) (Mo (09 erg) (Mo) (Mo) (10% erg/s (days 
Text = 1 X 109cm 2x 10? cm 1 x 10? cm 2 x 10? cm 

Standard resolution models 

T20 0.29 14.9 >1.4(7.0) >1.6 (6.1) »48(-25) »3.7(-2.0 20.71 0.56 14.5 (7.7) 12.0 (26.0 

ADO9x1 0.87 13.2 0.08 (1.6 0.12 (1.5) 0.57 (—0.21) 0.53 (—0.19 0.04 0.01 0.50 ( 0.28 3.3 (13.1 

AD20x1 0.92 20.9 0.22 (6.0) — 0.25 (5.6 1.8 (—1.4) 1.3 (C14 0.16 0.00 2.6 (0.83) 4.4 (27.8 

AD20x2 2.33 152 0.96 (7.9 1.3 (6.6 3.5 (—1.7) 3.1 (—1.6 0.15 0.15 4.1(2.0) 10.9 (26.8 

Lower resolution models 

T201 0.29 13.0 > 0.62 (8.2) >1.1 (7.7) »18(-32) »1.5(-2)9 > 0.25 m - = 

ADO9x11 0.87 14.4 0.11 (1.5) 0.21 (1.4) 0.57 (—0.20) 0.54 (—0.18 0.04 E - - 

BHdisk 0.0 12.2 > 0.088 (— > 0.3 (— > 0.088 > 0.037 i - 


Values with “>” denote that they are still increasing at the end of the simulation. 


accretion shock, which is formed along the disk surface. 
In the late phase, we also define the mass accretion rate 
to the black hole by 


Magn = 


T—TAH 


V —gpu* dsp. (13) 


'The definitions of the unbound matter and explosion 
energy are the same as those in Fujibayashi et al. (2021). 
We first define the specific binding energy and binding 
energy flux density of the matter as 


T, 

€bind ‘= TR a (1 + Emm) (14) 
pu 

fiina = -Té — pu’ (1 + Emin), (15) 


where T;* and T;? are the time-time and time-space com- 
ponents of the energy-momentum tensor and £j is the 
minimum specific internal energy of the employed EOS 
table (z —0.0013c? for SFHo and DD2). We then define 
the condition for unbound matter as e€ping > 0. We find 
however that the binding energy of the outer layer of 
the star in our computational domain is not accurately 
calculated due to a numerical error accumulated in long- 
term simulations. Therefore, we decide to calculate the 
diagnostic explosion energy as well as the mass of the 


unbound matter as 


— t 3 
Eexp = C€bind gu vy — gd x 
€bind 20,T feet 


i 
F J / fina dsk dt, (16) 
€bind >0, T =T ext 
Mej = pu /—gd? x 
€bind 2 0,T <Text 


t 
+ J / using dsp dt, (17) 
€pind 20,r—rext 


where ds; is the surface element of the sphere with the 
radius req. These quantities are defined by the volume 
integral for the matter of the positive binding energy 
inside an extraction radius rex, plus the time integral 
for the components of the positive binding energy flux 
at the radius. 

A part of the stellar matter is located outside the ex- 
traction radius and even outside the computational do- 
main. Its mass and binding energy can contribute to 
the ejecta mass and explosion energy. To estimate their 
contribution, we first obtain the enclosed mass at the 
extraction radius when the shock wave reaches the ex- 
traction radius. We then estimate the binding energy 
of the matter above the radius with the same enclosed 
mass in the pre-collapse profile. We compare the ex- 
plosion energies with different extraction radii in § 3.4. 


We note that the matter in the progenitor stars is 
composed mainly of He, !?C, 160, and Ne. When 
160 burns into ?9Ni, the rest-mass energy of 0.67 MeV 
per baryon is released into the internal energy, which 
can be an important energy source of the explosion of 
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the star (if the internal energy is not carried away by 
the neutrino emission). However, this effect is absent 
in our simulation because of the assumption of nuclear 
statistical equilibrium (NSE) in constructing the EOSs, 
for which the low-temperature matter in the computa- 
tional domain is assumed to be composed mostly of °° Ni. 
In addition, °°Ni is photo-dessociated into lighter nuclei 
if the temperature of the matter increases to Z 7 GK, 
consuming more internal energy than in the dissocia- 
tion of !?C or 160. This photo-dissociation spuriously 
takes place when a shock is formed near the disk around 
a black hole and propagates in the outer region of the 
star. T'hus, the explosion energy in our present simula- 
tions is underestimated. This point will be discussed in 
83.7. 


2.4. Tracer-particle method 


'To perform nucleosynthesis calculation, we apply our 
post-process tracer-particle method for the results of 
our simulations. The method is the same as that in 
Fujibayashi et al. (2022). The tracer particles are dis- 
tributed for 128 polar angles in the range of 0 = [0 : 1/2] 
on the arc with the radius of rex, = 2 x 10^ km. The 
particles are continuously set with the time interval of 
Atset := Tea4A0/(v"), where A0 = 7 /2/128 and (v") 
is the average radial velocity of the ejecta at the ex- 
traction radius. The mass of each particle is deter- 
mined based on the mass flux at the extraction radius 
as Am = rea? AQpu" /—gAtset, where AQ is the solid 
angle element. 


3. RESULT 
3.1. Evolution after bounce to disk formation 


For all the collapse models, a proto-neutron star is 
initially formed after the core bounce together with the 
standing accretion shock formation. Subsequently, due 
to the matter accretion from the outer region, the proto- 
neutron star collapses to a black hole. The black-hole 
formation time is tpb = 0.3, 2.3, 0.9, and 0.9 s for models 
T20, AD20x2, ADO9x1, and AD20x1, respectively (denoted 
by tsu in Table 2). Here tpp denotes the time after the 
core bounce. For model AD20x2, the formation of the 
black hole is delayed due to the significant centrifugal- 
effect associated with the rapid rotation. The earlier 
collapse for model T20 reflects a high mass-infall rate 
after the core bounce (see Fig. 2). These results illus- 
trate that the formation process of the black hole (and 
subsequent evolution process of the black hole and disk) 
depends on the profiles of the density and specific angu- 
lar momentum of the progenitor stars. 

For seconds after the black-hole formation, the in- 
falling matter does not have enough angular momen- 
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Figure 2. Top: Mass infall rate across the surface of r = rin 
(solid curves) and mass accretion rate onto the black hole 
(dashed curves). Middle: Total neutrino luminosity (solid 
curves) and total viscous heating rate (dashed curves) inside 
the shock wave. Bottom: Neutrino cooling efficiency defined 
by the total neutrino luminosity divided by the mass accre- 
tion rate onto the black hole. For model BHdisk, tpb = 0 
corresponds to the beginning of the simulation. 


tum to form disks, and thus, it is swallowed simply by 
the black hole. A geometrically thin disk starts forming 
at typ © 1.5, 5.0, 7.0, and 10.5s at which the black- 
hole mass is Afen AA ~ 3.5, 4.5, 5.0, and 7.0 for mod- 
els T20, AD20x2, ADO9x1, and AD20x1, respectively (see, 
e.g., panel (a) of Fig. 3 for model T20). The above black- 
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Figure 3. Snapshots at tpb» = 2.9 (a), 3.0 (b), 15.3 (c), and 23.4s (d) for model T20. The panels (a)-(d) show the snapshots 
just prior to the formation of a geometrically thick disk, in the quasi-steadily growing phase of the disk, at the onset of the 
outflow, and the final snapshot, respectively. Each panel has four sub-panels of rest-mass density (top-left), entropy per baryon 
(top-right), temperature (bottom-left), and electron fraction (bottom-right). The black solid circle for the first panel shows 
the region inside the apparent horizon (for other panels, it is too small to be seen because the plotted region is much wider 
than GMpy/c’). Note that the regions of the plots are different for each snapshot. An animation for this model is available at 
https://www2.yukawa.kyoto-u.ac.jp/~sho.fujibayashi/share/anim_20_hres_visc_multiscale.mp4. 


hole mass is consistent with the ones inferred from the 
distributions of density and specific angular momentum 
of the progenitor stars (see § 2.2). 

The disks at their formation is geometrically thin be- 
cause of the lower pressure inside the disk than the ram 
pressure of the infalling matter (Sekiguchi & Shibata 
2011). The disks then become geometrically thick when 
the pressure in the outer part of the disk and ram pres- 
sures of the infalling matter become comparable. The 
vertical expansion of the disk occurs at tpb 7 3.0, 10, 13, 
and 21s for models T20, AD20x2, ADO9x1, and AD20x1, 
respectively (see panel (b) of Fig. 3 and panels (a), (d), 
and (g) of Fig. 5). In the presence of the viscosity, the 


disk expansion occurs at time slightly later than that 
in the corresponding non-viscous model. The reason for 
this is that the viscous angular momentum transport ac- 
celerates the accretion of the disk matter onto the black 
hole and the increase of the pressure in the disk is de- 
layed. 

After the vertical expansion of the disk, a shock sur- 
face between the disk and infalling matter expands. This 
enhances the dissipation of the kinetic energy of the in- 
falling matter at the shock more efficiently. As a con- 
sequence, the neutrino luminosity increases during this 
phase. Indeed, the neutrino luminosity depends on the 
mass infall rate (see Fig. 2), and for rapidly rotating 
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Figure 4. Density (top-left), neutrino cooling rate density (bottom-left), cooling efficiency (top-right), and viscous heating 
rate density (bottom-right) for model T20 at tpp © 11.6 (left panel) and 15.4s (right panel), respectively. 


models T20 and AD20x2, it is higher because the mass 
accretion rate is higher. 

After the onset of the vertical disk expansion, rapidly 
rotating models T20 and AD20x2, and other models 
ADO9x1 and AD20x1 show different evolution process in 
terms of mass infall rate and neutrino cooling efficiency. 
We describe the evolution processes in the following sub- 
sections seperately. 


3.2. Models T20 and AD20z2: Models of higher 
infalling rate at disk formation 


For model T20, the disk evolves in a quasi-steady man- 
ner for ~ 10s after the onset of the vertical expan- 
sion. In this phase, the disk mass increases with a rate 
Maus = Miti — Man 2 0.2Mo /s (see top panel of Fig. 2 
for 4s S$ tpb S; 15s). This rate is much higher than other 
models, because for model T20, the disk starts forming 
at a relatively earlier stage of the collapse. 

In the quasi-steady phase of the disk, the neutrino 
cooling efficiency, defined by the neutrino luminosity L, 
divided by Mpyc?, is several to ten per cent (see the 
bottom panel of Fig. 2), which is sufficiently high for 
neutrinos to carry away the energy generated by viscous 
heating. That is, a neutrino-dominated accretion flow 
(NDAF) is established in this phase. The middle panel 
of Fig. 2 shows that the neutrino luminosity is ~ 10% 
higher than the viscous heating rate. This situation is 
possible because of the presence of an additional heating 
by the infalling matter, which releases its gravitational 
binding energy at shocks, generated on the disk surface, 


with the rate of 


L GMpuMaisk 
intal Y S c 
T'disk 
Mpu Maisk Tak Y ` 
x 1 51 
dale EH (3 LA i 


(18) 


where the infalling matter is assumed to settle into a 
Keplerian motion at a radius fqisk- Interestingly, the 
neutrino luminosity is approximately constant in time 
for the quasi-steady phase, 4s S tpb S 14s. 

After the infalling rate Mian starts decreasing at tpb Z 
14s, the viscosity-driven outflow is launched because the 
viscous heating plus shock heating dominate over the 
neutrino cooling on the surface region of the disk. Fig- 
ure 4 shows the rest-mass density, local neutrino cooling 
efficiency, neutrino cooling rate, and viscous heating rate 
before (tpb = 11.6s: left) and after (tpp = 15.4s: right) 
the onset of the outflow. Here, the local neutrino cool- 
ing efficiency is defined by q7 die, + q} ), where qz is 
neutrino cooling and heating rate densities and qi, is 
the viscous heating rate density (see Fujibayashi et al. 
2018 for its definition). The balance between the neu- 
trino cooling and the viscous heating is established in 
the mid-plane region (white regions) at both phases. On 
the surface of the disk, where the rest-mass density is 
~ 108 g/cm, a velocity shear between the disk matter 
and the infalling matter is clearly present and produces 
a high viscous heating rate (see the bottom-right plot 
in each panel). In addition, it is found that the neu- 
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Figure 5. The same as Fig. 3 but for models AD20x2 (top panels), ADO9x1 (middle panels), and AD20x1 (bottom panels). For 
each model, the panels from left to right show the snapshots at the formation of a geometrically thick disk, after the onset of 


the outflow, and the final snapshots, respectively. 


trino cooling is not efficient on the disk surface because 
of the low temperature there. This results in a launch 
of the outflow from the surface of the disk once the ram 
pressure of the infalling matter becomes weak. 

Model AD20x2 also exhibits a phase in which the vis- 
cous heating and neutrino cooling rates are comparable 
for tpb ~ 10-16s. Thus, the evolution process of the 
disk is similar to that of model T20; for model AD20x2, 
the explosion also sets in at ~ 10s after the formation 
of a geometrically thick disk. 


For rapidly rotating models T20 and AD20x2, a high- 
mass infall phase continues prior to the disk outflow for 
a long timescale of ~ 10s. The reason for this is that 
during such a phase, the ram pressure from the infalling 
matter is strong and the temperature in the shocked re- 
gion is high enough for enhancing the neutrino cooling. 
As a result, the thermal pressure generated by the vis- 
cous heating cannot be high enough for launching the 
disk outflow. 
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Figure 6. Mass histograms as a function of the electron 
fraction (top), entropy per baryon (middle), and expansion 
timescale (r/v^; bottom) at T = 5 GK of the tracer particles 
that experience temperature higher than 5 GK. 


After the onset of the outflow from the disk, an ex- 
panding shock with a slightly oblate shape is formed (see 
panel (c) of Fig. 3 for model T20 and panel (b) of Fig. 5 
for model AD20x2). During the shock expansion, the 
matter infall onto the black hole and disk still continues 
in the polar and equatorial directions because the out- 
flow from the disk is launched along the geometrically- 
thick-disk surface. The outflow towards the polar and 
equatorial directions is suppressed in different reasons. 
In the equatorial direction, the outflow is prohibited due 


to the presence of the disk. Near the polar axis, the in- 
falling matter that passed through the expanding shock 
surface converges toward the polar region. As a result, 
the ram pressure near the polar region is enhanced and 
becomes larger than that in the other direction. This 
prevents the outflow toward the polar direction. The 
final snapshots for these models also show that the out- 
flow from the disk expands in a diagonal direction (see 
panel (d) of Fig. 3 and panel (c) of Fig. 5; the matter 
with higher entropy, s/kg = 30, is the outflow compo- 
nent launched from the disk). 

'The top panel of Fig. 6 shows the mass histogram as a 
function of the electron fraction of the ejecta. Although 
a neutron-rich region with Y, < 0.2 is present in the 
inner mid-plane region of the disk reflecting the high 
density there (see Figs. 3 and 5), the lowest values of 
Y, of the ejecta is ~ 0.46 and 0.40 for models T20 and 
AD20x2, respectively. The reason for these high values 
is that the outflow is launched from the disk surface 
region, in which the electron degeneracy is not as high 
as that in the disk mid-plane and the value of Y, is close 
to 0.5. For model AD20x2, the lowest value of Y, for the 
ejecta is smaller than for model T20. The reason for this 
is that not only the matter on the disk surface but also 
somewhat inside the disk surface is ejected for model 
AD20x2. For this model, the neutrino cooling is always 
less efficient than that for T20. This leads to the ejection 
of the matter deeper inside the disk with lower electron 
fraction due to the viscous heating. This indicates that 
the electron fraction of the ejecta depends on inherent 
profiles of the density and specific angular momentum 
of the progenitor stars. 

The electron fraction for some of the ejecta can be 
higher than 0.5. This reflects the fact that the positron 
capture proceeds in a shorter timescale than the electron 
capture under the condition of low electron degeneracy 
and mildly high temperature kgT > 1 MeV. In this case, 
the positron capture n + et — p+ Ve is energetically 
more preferred than the electron capture because the 
mass difference between a free neutron plus an electron 
and a free proton (Mn + Me — Mp)c? is not negligible 
compared with the kinetic energy of electrons ~ kgT 
(see Just et al. 2022c and Beloborodov 2003). 

From the ejecta with Y, 2 0.5 and with T = 5GK, 
a substantial amount of ?9Ni can be synthesized. We 
discuss this topic in 83.7. 

'The middle panel of Fig. 6 shows the mass histogram 
as a function of entropy per baryon, s/kp. Typical val- 
ues of s/kp are found to be 10-20 for rapidly rotating 
models T20 and AD20x2, which are lower than for mod- 
erately rotating models (ADO9x1 and AD20x1). This is 
because the outflow is developed well (~ 10s) after the 
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formation of quasi-steady disks in a milder manner for 
the rapidly rotating models than for the moderately ro- 
tating models. This trend is also illustrated in mass his- 
togram of the expansion timescale (see the bottom panel 
of Fig. 6). The expansion timescales for the rapidly ro- 
tating models are typically Z 0.1s, which is longer than 
for the moderately rotating models. The entropy and 
expansion-timescale distributions for the rapidly rotat- 
ing models are similar to those for model BHdisk. 

For a rapidly rotating model T20, a disk with mass of 
~ 2Mo still remains at the termination of our simula- 
tions and the bulk of the disk is still cooled by neutrino 
emission efficiently. When the cooling efficiency drops, 
a stronger outflow may be launched (e.g., Fernández & 
Metzger 2013; Just et al. 2015; Fujibayashi et al. 2020a; 
Just et al. 2022a). Suppose that 1096 of the disk mat- 
ter becomes ejecta eventually, the matter with mass of 
O(0.1) Mg becomes the ejecta in a later phase. Since 
the outflow matter experiences a condition with tem- 
perature higher than 5 GK, a non-negligible amount of 
56Ni may be synthesized in the later-time ejecta of this 
model. 

On the other hand, the disk mass is ~ 0.2Mo for 
model AD20x2 and the cooling efficiency is low at the 
termination of the simulation. The electron fraction for 
the bulk of the disk matter is frozen to be Y, = 0.4-0.5, 
because of a long weak interaction (electron/positron 
capture) timescales. Therefore, the electron fraction of 
the matter expected to be ejected in a longer timescale 
is likely to be 0.4-0.5 (see also 8 3.6). 


3.3. Models ADO9z1 and AD20z1: Models of lower 
infalling rate at disk formation 


For models ADO9x1 and AD20x1, the viscous heating 
dominates over the neutrino cooling in the entire phases 
after the disk expansion (see the middle panel of Fig. 2). 
The neutrino cooling efficiencies for these models are 
also lower than those for models T20 and AD20x2 (see 
the bottom panel of Fig. 2). This stems from small mass 
infall rates onto the disk and black hole and from result- 
ing weaker ram pressure of the infalling matter after the 
formation of the disk for models ADO9x1 and AD20x1. 
The dominance of the viscous heating leads to an early 
launch of the outflow (several hundreds of milliseconds) 
after the disk expansion for these models. In contrast to 
models T20 and AD20x2, the disk mass and density do 
not increase significantly for models ADO9x1 and AD20x1 
before launching the outflow. As a result, the outflow in 
the equatorial direction is not suppressed by the disk but 
rather it is enhanced (see panel (e) and (h) of Fig. 5). 
The final snapshots (panels (f) and (i) of Fig. 5) also 


show the deformed profile of shock surfaces and the out- 
flow from the disk toward the equatorial direction. 

For these models, the value of Y, is not low (at lowest 
0.47; see the top panel of Fig. 6) because of their low disk 
density and low neutrino cooling efficiency throughout 
the disk evolution. This leads to weak electron degen- 
eracy of the disk matter and thus results in the higher 
value of Y,. Some components have the electron frac- 
tion even higher than 0.5 due to the same reason as for 
models T20 and AD20x2 (see § 3.2). The result of these 
models are similar to those found in Just et al. (2022a). 
Because of the high electron fraction of the ejecta, we 
do not expect an r-process nucleosynthesis (see § 4.2). 

'These models also have larger values of entropy per 
baryon s/kp ~% 30-50 and shorter expansion timescales 
r/v" ~ 50ms than rapidly rotating models T20 and 
AD20x2 (see the middle and bottom panels of Fig. 6). 
In our interpretation, this result stems from the follow- 
ing three reasons. First, for moderately rotating models, 
the typical radius of the disk is smaller at the onset of 
the outflow because of the shorter time of the onset of 
the outflow after the disk formation. This requires a 
larger heating efficiency for the disk matter to be un- 
bound and for the outflow matter to have the higher 
velocity. Second, since the outflow sets in before the 
disk settles into a quasi-steady phase, the shear of the 
poloidal velocity field in addition to that of Keplerian 
motion leads to more efficient viscous heating. Third, 
the entropy generated by the viscous heating is not effi- 
ciently lost by the neutrino emission because of the lower 
neutrino cooling efficiency. These properties of moder- 
ately rotating models produce the ejecta with a higher 
entropy and a shorter timescale. 

As these results illustrate, the evolution of the disk, 
the timing for the onset of the explosion, and the elec- 
tron fraction of the ejecta depend strongly on the dis- 
tribution of the specific angular momentum of the pro- 
genitors. For the models with relatively low angular 
momentum (ADO9x1 and AD20x1), the explosion occurs 
in a short timescale after the formation of a geometri- 
cally thick disk, at which the ram pressure of the in- 
falling matter is low enough for launching an outflow. 
In this case, the electron fraction of the ejecta result- 
ing from the explosion cannot be very low. By contrast, 
for the models with relatively high angular momentum 
(T20 and AD20x2), the explosion takes place at ~ 10s 
after the formation of the disk, and in this case, the elec- 
tron fraction of the ejecta can be low with Y, < 0.45, 
and thus, light trans-iron nuclei can be synthesized (see 
8 4.2). 


3.4. Ejecta mass and explosion energy 
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Figure 7. Time evolution of ejecta mass (top) and ex- 


plosion energy (bottom) as functions of post-explosion time 
t — texp. The dashed and solid curves are the values calcu- 
lated with rext = 1 x 10? cm and 2 x 10? cm, respectively. 
The values in the legend denote the mass (in units of solar 
mass) and binding energy (in units of 10?! erg) above the ex- 
traction radius when the shock wave reaches these radii (see 
$2.3). 


Figure 7 shows the mass of the ejecta and a diagnostic 
explosion energy as functions of time after the explosion, 
t—texp. Here, the time of the explosion is defined as the 
time at which the explosion energy reaches 1 x 10°° erg. 

For all the viscous-hydrodynamics simulations, we find 
the explosion by the energy injection from the disk 
around the black hole. For rapidly rotating progenitor 
models T20 and AD20x2, after the launch of the out- 
flow, the explosion energy increases beyond 1 x 10?! erg 
with a timescale of ~ 1s. For these models, the ex- 
plosion energy at the termination of the simulation 
reaches > 3 x 105! erg, and is still increasing. Since 
the timescale of the energy increase appears long, we 
stopped the simulation by t = 30s. Thus, the ejecta 
mass and explosion energy which we describe below give 
the lower limit. 

At the termination of the simulations, disks with ~ 2 
and 0.2Mo are present around the black holes for mod- 
els T20 and AD20x2, respectively. In addition, matter in- 


fall still continues around the central region. For these 
models, therefore, an energy budget to provide more 
explosion energy is still present. These results suggest 
that if the massive progenitors are rapidly rotating, a 
high-energy supernova-like explosion could occur from 
the disk outflow triggered by viscous heating. We dis- 
cuss this point in 84.1. 

For models ADO9x1 and AD20x1, for which the angu- 
lar momentum of the progenitor stars is relatively low, 
the formation of the disk is delayed, and the onset of 
the explosion occurs in a late phase (see Table 2). As a 
consequence, only a small amount of mass remains out- 
side the black hole at the explosion. The small energy 
budget at the formation of the disk leads to a slow in- 
crease in the explosion energy. For these models, the 
explosion energies amount to œ~ (0.5-1) x 10?! erg at the 
termination of the simulations, which are comparable 
or smaller than the canonical value for core-collapse su- 
pernovae. The ejecta masses for these models are 0.1— 
0.3Mo at the termination of the simulations, which are 
an order of magnitude smaller than that for canonical 
supernovae. 

For all the collapse models, appreciable stellar matter 
is still present above the extraction radius. This can 
contribute additionally to the ejecta mass and explosion 
energy. The masses outside the extraction radius are 
larger than those of the ejecta, and thus, they may have 
an impact on the ejecta velocity. The terminal aver- 
age velocity of the ejecta depends on the amount of the 
stellar matter that becomes ejecta in the later phase. 

The mass and binding energy of the matter outside 
the extraction radii rex, = 1 x 10? cm and 2 x 10? cm 
are listed in Table 2 (see the values in parentheses for 
columns of Mej and Exp; also see Fig. 7). At the ter- 
mination of the simulations, the absolute values of the 
binding energy are below the explosion energies for mod- 
els T20, ADO9x1, and AD20x2. Thus, the explosion is 
likely successful for these models. For model AD20x1, 
on the other hand, the value is slightly above the explo- 
sion energy for Text = 2 x 10? cm. Thus, it is not clear 
whether this model results in successful explosion. Note 
however that the bound matter outside the extraction 
radius at the termination of the simulation has a suf- 
ficiently large angular momentum to circularize around 
the black hole (see Fig. 1). Thus, a possible increase 
may be expected in the explosion energy for the later 
phase if such a matter falls into the central region to 
power an additional outflow. We are planning to per- 
form a longer-time simulation in our future work to ad- 
dress the possibility of the explosion for these models. 
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As already mentioned in 8 2.3, the composition of the 
computational domain is assumed to be that in NSE in 
our simulation. In reality, the outer region of the star 
is composed mainly of !?C, 160, and ??Ne, and hence, 
the nuclear burning of them in the stellar mantle swept 
by the shock wave can provide an additional energy. In 
8 3.7, we will investigate this more quantitatively. 


3.5. Resolution study 


In Fig. 8, the results for lower resolution models T201 
and ADO9x11 are compared with their corresponding 
standard resolution models. The results for models 
T201 and ADO9x11 agree qualitatively with those of the 
higher-resolution models T20 and ADO9x1. However, 
a quantitative difference is found in the onset time of 
the outflow, which is tpp ~ 12.3s and 14.5s for T201 
and ADO9x11, while tpb ~ 15.0s and 13.3s for T20 and 
ADO9x1. The fluid dynamics after the onset of the out- 
flow is also somewhat different, and thus, the time evo- 
lution of the ejecta quantities is quantitatively different. 
Despite these differences, however, the saturated values 
of ejecta mass, its Y, distribution, and explosion energy 
for model ADO9x11 agrees well with those for ADO9x1. 
Also for model T201, the Y, distribution agrees well with 
model T20. For the lower-resolution run (T201), the ex- 
plosion energy is lower than the higher-resolution run 
(T20). This suggests that the explosion energy might be 
underestimated in our present studies. 


3.6. Comparison of the collapse models with model 
BHdisk 


For model BHdisk, which consists of a 3Mo disk 
around a 10M black hole with the dimensionless spin 
of 0.6, the disk matter accretes onto the black hole in 
a quasi-steady manner in the first 5 seconds. The out- 
flow is then launched at t ~ 5s because the neutrino 
cooling rate has dropped far below the viscous heating 
rate (see the middle panel of Fig. 2). The increase rate 
in the explosion energy for model BHdisk is much lower 
than for the collapsing star models (see bottom panel 
of Fig. 7). The primarily reason for this is that, for 
the collapsing stars, a velocity shear is present not only 
in the accretion disk with the nearly Keplerian motion 
but also on the surface of the disk resulting from the 
infalling matter onto the disk, which is absent for model 
BHdisk. The strong shear on the disk surface signifi- 
cantly enhances the viscous heating rate, which results 
in the higher increase rate of the explosion energy. 

The origin of the velocity shear associated with the 
infalling matter is different from the Keplerian motion 
of the disk. Here, for the latter we suppose that the MRI 
turbulence (Balbus & Hawley 1991) is the origin of the 


effective viscosity. For the surface region of the disk, 
on the other hand, we suppose that the shear region 
on the disk surface should induce the Kelvin-Helmholtz 
instability. In such regions, magnetic fields are supposed 
to be enhanced significantly leading to the development 
of a turbulence and dissipating the kinetic energy of the 
infalling matter (e.g., Zhang et al. 2009; Obergaulinger 
et al. 2010; Viganó et al. 2020). Thus, it is natural 
to consider that the effective viscosity on such a region 
is also high.? However, to clarify this process in the 
first-principle way, magnetohydrodynamics simulation is 
necessary. In the future work, we plan to perform this 
to confirm that our assumption is indeed correct. 

For model BHdisk, we stopped the simulation at t ~ 
21s. For this model, the mass ejection still continues 
with the ejection rate slightly higher than the mass ac- 
cretion rate onto the black hole. This suggests that 
a large fraction of the disk matter will be eventually 
ejected from the system. At the termination of the sim- 
ulation, the total ejecta mass is only ~ 0.1 Mo. However, 
the disk mass is still zz 2M. Extrapolating the mass 
ejection rate at the final time ~ 10^? M/s, we infer 
that the mass ejection continues for more than 100s for 
this model. 

'The mass ejection for model BHdisk sets in after the 
neutrino cooling efficiency of the bulk of the disk drops 
(cf. the middle panel of Fig. 2). The ejecta for model 
BHdisk has a low-Y, component down to Y, z 0.4. 
'The value is determined by the electron fraction in elec- 
tron/positron capture equilibrium when the timescale of 
these reactions becomes comparable to that of the vis- 
cous expansion (see Fujibayashi et al. 2020c; Just et al. 
2022a). Thus, the mass ejection mechanism for rapidly 
rotating models is qualitatively similar to that for model 
BHdisk. 


3.7. Production of ?9 Ni 


Using the time evolution of the temperature and den- 
sity along tracer particles, the post-process nucleosyn- 
thesis calculations are performed with the nuclear reac- 
tion network code "NET (Wanajo et al. 2018). The initial 
composition of the nucleosynthesis calculation depends 
on the thermal history of the tracer particles. If the 
maximum temperature along a particle is higher than 
10GK, we start the nucleosynthesis calculation at the 
time that the temperature decreases to T' 2 10 GK with 


? For this case, the kinetic viscosity is likely to be proportional to 
the infall velocity, not to the sound velocity, as v = fturVintall; 
where Gaata is comparable to or larger than cs, and thus, our 
present treatment for the viscous coefficient may be conservative 
if &ur/H = O(10?). 
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for models T20, T201, ADO9x1, and ADO9x11. 


the mass fraction of free protons and nucleons, Y, and 
1 — Y, respectively. For the tracer particles with the 
maximum temperature lower than 10 GK, we start the 
nucleosynthesis calculation at t = 0 of the simulation 
with the composition depending on the position of the 
particle in the progenitor star (mostly consisting of 160 
and ?°Ne). 

The resulting mass of ?9Ni is listed in Table 2. It 
is found that for models T20 and AD20x2 the mass of 
56Ni amounts to ~ 0.56 and 0.15Mo, respectively, re- 
flecting the large mass of the ejecta that experiences 
high temperature > 5 GK. The amount of ?9Ni for the 


model T20 is high enough for explaining a typical ?9Ni 
mass inferred from the observations of type Ic-BL super- 
novae z 0.4Mo (see, e.g., Cano et al. 2017 and Anderson 
2019). It was pointed out that the ?9Ni mass inferred 
with the so-called *Arnett's rule" (Arnett 1982) is pos- 
sibly overestimated (e.g., Meza & Anderson 2020, who 
suggest that the °°Ni masses inferred from the radioac- 
tive tail luminosity for nearby two type Ic-BL super- 
novae, SN2009bb and SN2016coi, are 0.08 and 0.10Mo, 
respectively). Thus, the mass of °°Ni found for model 
AD20x2 could still be high enough for explaining those 
supernovae. These results suggest that very massive and 
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rapidly rotating stars leading to black hole and massive 
disk formation are candidates for the progenitors of type 
Ic-BL supernovae. 

We also note that the numerical simulations for mod- 
els T20 and AD20x2 underestimate the total ejecta mass 
because we stopped the simulations at a time when the 
ejecta mass and explosion energy are still increasing. 
Our results here indicate the lower bound for the ?9Ni 
mass. 

For models ADO9x1 and AD20x1, the °°Ni masses are 
smaller than 0.1 Mg (0.01 and 0.06Mo, respectively), 
because of their smaller ejecta masses. Thus for moder- 
ately rotating progenitor models, the ?9Ni mass is likely 
to be comparable to or smaller than that for an ordinary 
supernova. These models predict the presence of mod- 
erately bright, but rapidly varying optical transients as 
found in 8 4.1. 

It is important to note that, as can be found in Table 2, 
the produced °°Ni mass fraction relative to the mass of 
ejecta exceeding 5 GK varies from 2496 (AD20x2) to 7896 
(T20) depending on the electron fraction, entropy, and 
expansion timescale of the outflowing matter for each 
model. The conditions of Y, = 0.5, low entropy, and 
slow expansion are favored for the efficient production 
of ?9Ni (T20 is the case; see Fig. 6). Therefore, the mass 
of ejecta with > 5 GK (as frequently used in the litera- 
ture) only serves as a loose upper limit for the produced 
amount of ?9Ni. 

Suppose that ?9Ni are synthesized from !9O, the rest- 
mass energy released into the internal energy due to 
the nuclear burning is 7.0 x 1099, 1.6 x 104°, 7.5 x 1049, 
and 1.9 x 109? erg for models T20, ADO9x1, AD20x1, and 
AD20x2, respectively. These contributions are several 
to tens percents of the explosion energy estimated in 
8 2.3 and can have notable effects. Especially for model 
AD20x1, the explosion energy plus binding energy above 
the extraction radius rex, = 2 x 10? cm becomes nega- 
tive but the value is comparable to the energy gener- 
ated by nuclear burning. Thus, to clarify whether such 
a marginal model explodes successfully, feedback of the 
nuclear reaction has to be taken into account in hydro- 
dynamics simulations (see Bollig et al. 2020 for a recent 
attempt). 


4. DISCUSSION 
4.1. Optical transients 


Using the ejecta properties obtained in the present 
study, we analytically calculate bolometric luminosity 
models for photons following Arnett (1982). The ther- 
malization efficiency for gamma-rays is estimated follow- 
ing Colgate et al. (1997) with the optical depth for non- 
thermal gamma-rays ky = 0.03cm?/g. For the optical 
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Figure 9. Bolometric light curves for standard reso- 


lution models. The time origin is chosen to be the peak 
time for each curve. The solid curves denote the light 
curves with ejecta mass and explosion energy extracted for 
Text = 2 x 10? cm (see Table 2 for the values). The dashed 
curves denote the ones with the assumption that a half of 
the mass and binding energy outside the extraction radius 
contribute to the ejecta properties. The point shows the time 
at which the ejecta becomes optically thin, T = &pR = 1, to 
thermal photons for each model. The shaded regions denote 
the templates of the bolometric light curves with standard 
deviations for type Ib, Ic, and Ic-BL supernovae taken from 
Lyman et al. (2016). 


depth of thermal photons, we simply set & — 0.1 cm?/g. 
We note that the ^Arnett" model tends to infer a larger 
96Ni mass by a factor of a few than that inferred by the 
"radioactive tail" luminosity of a supernova, the latter 
being less ambiguous (e.g., Meza & Anderson 2020, Af- 
sariardchi et al. 2021, and Rodríguez et al. 2022). This 
indicates that the luminosity predicted by the Arnett 
model for a given ?9Ni mass may be underestimated 
by a factor of a few (see Dessart et al. 2015, 2016 and 
Khatami & Kasen 2019). We also note that for rapidly 
rotating models (T20 and AD20x2), light trans-iron ele- 
ments could be synthesized in the ejecta (see 8 4.2), and 
hence, the opacity for optical wavelengths may be higher 
than 0.1 cm?/g. For more quantitative study, we obvi- 
ously need a radiation transfer simulation for photons 
taking into account a realistic opacity table. 

Figure 9 shows the bolometric light curves for all the 
models (with standard resolution) investigated in this 
paper. For deriving the solid curves of Fig. 9, we take 
into account only the ejecta mass and explosion energy 
extracted at re = 2 x 10? cm. For models ADO9x1 and 
AD20x1, the luminosity evolves rapidly. The rise times, 
defined by the time until the maximum luminosity is 
reached, are trise % 3.3 and 4.4d, respectively, for these 
models (see Table 2). Such fast transients may be dis- 
covered in the future high-cadence transient surveys. On 
the other hand, for models T20 and AD20x2, trise % 10d, 
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which is consistent with that of type Ib/c supernovae 
(see, e.g., Taddia et al. 2018). This indicates that there 
could be a possible subclass of type Ib/c supernovae 
driven by the disk outflow in the black hole-forming core 
collapses of rotating massive stars. 

For model T20, we note that the explosion energy, 
ejecta mass, and the mass of °°Ni do not saturate at 
the termination of the simulation. This indicates that 
the luminosity and timescale shown here are underesti- 
mated. Therefore, this model would represent more en- 
ergetic transients than typical supernovae such as type 
Ic-BL (see also explosion energy for the similar scenario 
found in Just et al. 2022b). 

'To investigate the possible effects of the stellar mat- 
ter outside the extraction radius, we also calculate the 
light curve for each model assuming that a half of the 
mass and binding energy of the matter for r > rex, = 
2 x 10? cm contribute to the ejecta mass and explosion 
energy (see the dashed curves in Fig. 9). Because of the 
increased ejecta mass and decreased explosion energy, 
peak luminosity and timescale for all the light curves 
become smaller and longer, respectively. 

For rapidly rotating models T20 and AD20x2, the 
timescales of the light curves become longer by a fac- 
tor of ~ 2 and are longer than those typically found for 
type Ib/c supernovae. Karamehmetoglu et al. (2022) 
reported such long-timescale type Ib/c supernovae re- 
cently. The long-timescale supernovae are reported as 
the explosions of massive (Mzams Z 25Mo) stars with 
the explosion energies comparable to typical type Ib/c 
supernovae, but they are infered to have a larger amount 
(> 0.1Mọ) of DN) mass. These facts indicate that 
our rapidly rotating models can explain such subclass 
of type Ib/c supernovae. 

For a moderately rotating models ADO9x1 and AD20x1, 
on the other hand, the timescales becomes ~ 13 and 
28d, respectively, comparable to or longer than that 
of typical type Ib/c supernovae. The peak luminosity 
is, however, about 10 times dimmer than that of typi- 
cal type Ib/c supernovae. This indicates that rotating 
massive stars exploded by outflows from black hole-disk 
system can produce a variety of transients depending on 
the rotation profile of the progenitors and the presence 
of the stellar envelope, although the explosion mecha- 
nism is qualitatively universal. 

As found in the comparison of the solid and dashed 
curves in Fig. 9, the features of the bolometric light 
curves depend on the possible contribution of the mat- 
ter outside the extraction radius. The quantitative pre- 
diction of the optical transients requires us to perform 
simulations for entire stars until the outer layer of the 
star is swept by the shock wave. 


Another possible astrophysical transients can be pow- 
ered by the interaction of the ejecta with a circum-stellar 
medium that can result from the strong mass loss of their 
progenitor prior to the stellar core collapse. The progen- 
itor models provided by Aguilera-Dena et al. (2018) are 
likely to be surrounded by a dense, massive (~ 0.1-1Mo 
within ~ 10! cm) circum-stellar medium at the core col- 
lapse. Since the ejecta mass is an order of 0.1-~ 1M, 
for models ADO9x1 and AD20x1, the ejecta will be signif- 
icantly decelerated in the circum-stellar medium, releas- 
ing a substantial fraction of its kinetic energy < 10?! erg. 
'The optical depth of the circum-stellar medium is esti- 
mated as 
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Here, the opacity for photons is assumed to be dom- 
inated by the Thomson scattering of fully ionized 
medium with Y, — 0.875 (ie., hydrogen and helium 
with mass fractions 0.75 and 0.25, respectively). Since 
the circum-stellar medium is optically thick, the released 


energy diffuses out from the circum-stellar medium with 
the diffusion time. The luminosity is then estimated as 


Lo E€Eexp . Resm 


taif ` Uoclaift 


=i 
x 5.5 x 1043 NE C 
RAUTEE (sas cm?/g 


2 
€ Uoo Resm 
2 
i (5) (; x dë ES SS 


where v4, = J/2E V / (Me + Mcsm) is the terminal ve- 
locity of the ejecta plus circum-stellar medium, e is the 
radiation efficiency, and 


"T + Mcsm) 
Laun nd EEN 


Jan 
1/2 
K 
zx 26d| ———5s— 
(oor) 
" Mej + Mosm Sr exp TNT (21) 
2M 105! erg f 


is the diffusion time of the expanding ejecta plus circum- 
stellar medium (see, e.g., Matsumoto & Metzger 2022 for 
similar expression). The second factor of the first line 
in Eq. (20) is the contribution of the adiabatic cooling. 
This can naturally lead to an optical transient like su- 
perluminous supernovae (see, e.g., Moriya et al. 2018 for 
a review). We note, however, that the properties of the 
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transient depend not only on the mass and radius but 
also on the density profile of the circum-stellar medium 
(see, e.g., Chevalier & Irwin 2011 and Suzuki et al. 2020 
for the curcum-stellar medium like a stationary wind). 


4.2. Possible synthesis of light trans-iron nuclei 


Figure 10 shows isobaric mass fractions obtained 
by the nucleosynthesis calculations for our standard- 
resolution models. We find prominent peaks at alpha- 
nuclei (with A multiple of 4) and A — 56. The peaks 
at A = 12 and 16 reflect the initial composition of !?C 
and !9O, respectively, in progenitor stars. The bulk of 
these nuclei remains unprocessed owing to relatively low 
temperature achieved. On the other hand, the nuclei at 
peaks of A — 20-40 (corresponding to elements Ne, Mg, 
Si, S, Ar, and Ca) are synthesized from !?C and 160 
in the ejecta that experience higher temperature but 
lower than that required for achieving NSE. The peak 
at A = 56 corresponds to ?9?Ni, which is synthesized 
predominantly in NSE with Y, = 0.5. Interestingly, a 
certain amount of nuclei heavier than the iron group 
(A = 60), up to A z 90, is found to be synthesized for 
a rapidly rotating model AD20x2. It is known that in 
slightly neutron rich (Y, = 0.4) ejecta, such trans-iron 
nuclei are synthesized predominantly in quasi-nuclear 
statistical equilibrium (QSE, Meyer et al. 1998; Wanajo 
et al. 2018) under an alpha-rich condition. Since the 
neutron-richness is not very high, heavy r-process ele- 
ments with A > 100 are not synthesized for any of the 
present models. 

'The first peak nuclei of r-process, especially Zr and Y 
(synthesized in QSE here), are known to have opacities 
higher than those of iron-group elements (Kawaguchi 
et al. 2021). Therefore, if such elements are appreciably 
synthesized, the resulting optical transients may have 
longer timescales than those without such elements. In 
addition, the peak luminosity will be lower and the spec- 
trum could be redder. To quantify the light curve and 
spectrum, we plan to perform a radiation transfer sim- 
ulation in follow-up work. 


4.3. Implications for gamma-ray bursts 


Figure 11 shows the masses (top) and dimensionless 
spins (bottom) of the black holes for our standard res- 
olution models as functions of post-bounce time. Here, 
the mass and dimensionless spin of the black holes are 
estimated from the equatorial and polar circumference 
radii of the apparent horizon (e.g., Shibata 2016). 

For models T20 and AD20x2, the mass accretion 
from the disk onto the black hole continues for a long 
timescale (Z, 20s) because of the presence of the long- 
term high-mass-accretion rate phase. This leads to 
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Figure 10. Isobaric mass fraction for standard resolution 
models. Nuclei at A — 12 and 16 are, respectively, predomi- 
nantly unprocessed !?C and !$O in the progenitor stars. 
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Figure 11. Masses (top) and dimensionless spins (bottom) 
of black holes for our standard resolution models as functions 
of post-bounce time. 


rapidly spinning black holes, the dimensionless spin of 
which is 0.96 and 0.80, respectively, at the termination 
of the simulations. During the long-term evolution pro- 
cess of the black holes, the accretion disk is likely to be 
in a turbulent state associated with magnetohydrody- 
namical instabitlities such as the MRI (Balbus & Hawley 
1991) in a realistic situation. By this, the magnetic-field 
strength should be enhanced, and the magnetic flux pen- 
etrating the black hole is increased as a result of the mass 
accretion (and magnetic flux accretion) onto the black 
hole. Such a highly spinning black hole penetrated by 
a strong magnetic field could be a promising central en- 
gine for relativistic jet via the Blandford-Znajek process 
(Blandford & Znajek 1977) that is believed to power 
gamma-ray bursts (e.g., Gottlieb et al. (2022a)). 
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For models ADO9x1 and AD20x1, in which we use the 
angular momentum distribution of the original stellar 
evolution simulations, the mass accretion and thus the 
spin up of the black hole are suppressed because of 
the late-time disk formation and the quicker launch of 
the outflow from the disk. The dimensionless spins 
for these models are ~ 0.4 and 0.7 respectively. The 
Poynting luminosity for these models may not be large 
enough for launching intense electromagnetic waves via 
the Blandford-Znajek mechanism because of its strong 
dependence on the black-hole spin. Moreover, because 
of the low density of the disk for these models, strong 
magnetic fields may not be sustained in the vicinity of 
the black hole. This indicates that more rapidly rotat- 
ing progenitors than those of ADO9 and AD20 would be 
preferred for generating powerful relativistic jets. 

It is known that at least a fraction of long-duration 
gamma-ray bursts are accompanied by type Ic-BL su- 
pernovae (e.g., Cano et al. 2017). As found in § 3.7, 
rapidly rotating models synthesize a large amount of 
°6Ni (> 0.1Mg), which is consistent with required 
amount to explain the light curves of type Ic-BL super- 
novae. In 8 3.4, it is also found that the rapidly rotat- 
ing models show a large average velocity of the ejecta 
> 2x10?cm/s, which is necessary for the broad-line 
features for these supernovae. In addition, the large ex- 
plosion energy > 105? erg observationally inferred from 
type Ic-BL supernovae is likely achieved if we consider 
longer evolution of the system than simulated in this 
study. Therefore, the rapidly rotating models in this 
study could reasonably represent the supernovae accom- 
panying long-duration gamma-ray bursts. 


5. SUMMARY 


In this paper, we studied the explosion in the rotating 
massive-star collapse leading to a black hole and a mas- 
sive disk in fully general relativistic radiation-viscous- 
hydrodynamics simulations with an approximate neu- 
trino radiation transfer. We employed evolved stars with 
compact core (Takahashi et al. 2018 and Aguilera-Dena 
et al. 2020) as the initial conditions for our numerical 
simulations. For one of the progenitor models, T20 (from 
Takahashi et al. 2018), a rapidly rotational profile is as- 
sumed, while for the other models ADO9 and AD20 (from 
Aguilera-Dena et al. 2020), the original or doubled angu- 
lar velocity is adopted. For all the models investigated in 
this paper, we found the formation of an accretion disk 
after a proto-neutron star collapsed into a black hole 
although the time at the onset of the disk formation de- 
pends strongly on the rotational profile of the progenitor 
stars. The evolution after the disk formation was also 


found to depend strongly on the degree of rotation of 
the progenitor stars. 

For models with rapidly rotating progenitors T20 and 
AD20x2, for which the mass infall rates to the disk and 
the black hole are high (Mian 2 0.3Mo/s), the disk 
evolves quasi-steadily prior to the onset of the outflow 
for ~ 10s after the formation of the disk. The outflow is 
launched from the surface of the disk after the mass in- 
fall rate decreases (i.e., the ram pressure of the infalling 
matter drops). At the launch of the outflow, the neu- 
trino cooling efficiency around the disk surface is lower 
than that deeper in the disk and the viscous heating 
rate. The ejecta mass and explosion energy amount to 
> 1Mg and = 3 x 10° erg, respectively. For model 
T20, they are still increasing at the termination of the 
simulations and the explosion energy may be as large 
as ~ LUTZ erg, The high explosion energy can explain 
the large kinetic energy inferred for type Ic-BL super- 
novae. The electron fraction of the ejecta is at lowest 
& 0.4 (for model AD20x2), which is not sufficient for an 
r-process nucleosynthesis. Indeed, the nucleosynthesis 
calculation showed that heavy nuclei are synthesized at 
most up to A ~ 100. However, there is low-electron frac- 
tion (Y; < 0.2) matter deep inside the disk in which the 
rest-mass density is high enough to enhance the electron 
degeneracy. If such a matter could be ejected by a very 
efficient mass ejection process with a shorter timescale, 
e.g., by magnetohydrodynamics processes, r-process el- 
ements might be synthesized. Magnetohydrodynamics 
simulation for this problem is one of our future issues. 

For model T20, the synthesized °°Ni mass (> 0.5Mo) 
as well as the predicted bolometric light curve is consis- 
tent with those observationaly inferred for type Ic-BL 
supernovae. On the other hand, the amount of ?9Ni for 
AD20x2 is ~ 0.15 Mo. The bolometric light curve for this 
model is suitable for a light-curve model of type Ib/c su- 
pernovae. This indicates that there might be a possible 
subclass of stripped-envelope supernovae driven by the 
outflow from a massive disk around a rapidly spinning 
black hole as a result of the collapse of a massive rotat- 
ing star. The possible existence of high-opacity trans- 
iron elements (such as Y and Zr) may lead to a longer 
timescale and redder transient. 

For moderately rotating models, ADO9x1 and AD20x1, 
for which the original angular momentum profiles ob- 
tained in stellar evolution calculations are employed, 
the rest-mass density and the neutrino cooling efficiency 
were already low at the disk formation. As a result, the 
outflow is launched before a massive disk is formed, and 
the explosion occurs at several hundreds of millisecond 
after the disk formation. Because the disk starts form- 
ing at a late stage of the stellar collapse (at tp, > 10s), 
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the mass of the envelop is relatively small, and hence, 
the energy budget is smaller than for the rapidly ro- 
tating models. As a consequence, the ejecta mass and 
explosion energy are relatively small as ~ 0.1Mọ and 
~ (0.5-1) x10?! erg, respectively. The electron fraction 
of the ejecta is always higher than 0.47, and thus, an 
r-process nucleosynthesis cannot proceed in the ejecta 
for these models. 

Depending on the effect of the mass in the outer layer 
of the star, the resulting optical transient can have a 
very short timescale (a few days) or that comparable 
to normal supernovae. If dense and massive circum- 
stellar media are present as predicted in Aguilera-Dena 
et al. (2018), a very bright (~ 104^ erg/s) transient with 
a timescale of months is expected due to the interaction 
of ejecta with the circum-stellar medium. 

To more rigorously predict observational features 
(photometric luminosity and spectra) of optical tran- 
sients based on our present results, we need to perform 
a detailed photon-radiation transfer simulation. The in- 
clusion of high-opacity trans-iron elements, which could 
exist in the ejecta, may drastically change the observa- 
tional feature, which will be investigated in our future 
work. 


We employed viscous hydrodynamics to take into ac- 
count angular momentum transport and the dissipation 
of kinetic energy to internal energy in the region in which 
a velocity shear or differential rotation is present. This 
enables us to approximately capture the effective vis- 
cosity induced by the turbulence due to magnetohydro- 
dynamical instabilities. However, it is obviously nec- 
essary to perform first-principle magnetohydrodynam- 
ics simulations for strictly exploring the effects of the 
angular momentum transport and (effectively) viscous 
dissipation. We plan to perform a three-dimensional 
radiation-magnetohydrodynamics simulation in our sub- 
sequent project. 
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